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BAYESIANS  EXTENSIONS  TO  A  BASIC  MODEL  OF  SOFTWARE  RELIABILITY 


by 

William  S.  Jewell 


1.  INTRODUCTION 

A  convenient  method  of  surveying  the  problems  of  concern  in  software 
reliability  is  to  scan  the  articles  in  the  August,  1979  issue  of  the  IEEE 
Transactions  on  Reliability  (Volume  R-28:3).  In  addition  to  various  relia¬ 
bility  management  and  computer  science  issues,  it  refers  to  the  process  of 
identifying  and  removing  computer  software  errors  through  program  testing. 

One  of  the  earliest  and  simplest  models  for  describing  the  stochastic 
process  of  error  occurrence  is  that  proposed  by  Jelinski  and  Moranda  [2]. 
Essentially,  they  assume  that  a  fixed,  but  unknown  number  of  defects  or  "bugs 
in  the  program  are  "competing"  equally  and  randomly  for  discovery,  each  at 
the  same  average,  unknown  rate.  The  objective  of  the  analysis  is  to  estimate 
the  number  of  undetected  faults  remaining  in  the  program  after  a  certain  time 
and/or  failure  history,  or  to  make  other  performance  guarantees  (Littlewood, 
[9]). 

Parenthetically,  we  should  note  that  many  variations  to  this  basic  model 
have  been  proposed — Littlewood  and  Verrall  [7]  [8],  Goel  and  Okumoto  [1],  and 
Langberg  and  Singpurwalla  [6]  analyze  those  competing  models  and  show  that 
they  can  be  cast  in  a  unified  framework.  More  recent  elaborations  are  in 
Moranda  [12],  Littlewood  [10],  and  Kreraer  [5]. 

In  this  paper  we  shall  stick  with  the  basic  Jelinski-Moranda  (J-M)  [2] 
model  for  simplicity,  and  to  highlight  our  Bayesian  extensions.  At  first, 
our  paper  is  tutorial,  following  the  Bayesian  development  of  Meinhold  and 
Singpurwalla  (M-S)  [11].  However,  the  following  important  extensions  are 


(a)  The  testing  protocol  is  permitted  to  run  for  a  fixed  length  of 
time — possibly,  but  not  necessarily,  coinciding  with  a  failure 
epoch; 

(b)  The  distribution  of  the  unknown  number  of  defects  is  generalized 
from  the  one-parameter  Poisson  distribution  by  assuming  that 
the  parameter  is  itself  a  random  quantity  with  a  Beta  prior 
distribution; 

(c)  Although  the  calculation  of  the  posterior  distributions  of  the 
parameters  leads  to  complex  expressions,  we  show  that  the  compu¬ 
tation  of  the  predictive  distribution  for  undetected  errors  is 
straightforward; 

(d)  Although  it  is  now  recognized  that  the  MLE's  for  reliability 
growth  can  be  very  unstable,  we  show  that,  if  a  point  estimator 
is  needed,  the  predictive  mode  is  easily  calculated  without 
obtaining  the  full  distribution  first. 

The  paper  concludes  with  a  numerical  example. 
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2.  BASIC  MODEL 

Suppose  there  are  N  defects  in  a  given  software  product  and  that,  as 
a  result  of  exercising  the  program,  failures  are  observed  to  occur  at  epochs 
t^j  <  t^)  <  ...  •  Each  failure  initiates  a  debugging  action  that  perma¬ 

nently  eliminates  the  error  causing  the  failure;  for  simplicity  we  assume 
that  the  time  to  find  and  fix  the  error  is  negligible,  or,  what  is  the  same 
thing,  we  measure  only  the  duration  of  time  while  running  the  program.  If 
ve  then  assume: 

(a)  The  overall  failure  rate  of  .the  program,  p(t)  hours  *  ,  between 

t^^  and  «  is  a  constant  equal  to  (N-i)0  hours  *  » 

(i  =  0,  1,  2,  ...  N-l  ;  t^j  =  0)  ,  where  0  is  a  unit  failure 
rate',  and, 

(b)  given  N  and  0  ,  the  inter- failure  times ,  ^  =  t^  -  t^+1j 

(i  =  0,  1,  ...  N-l)  ,  are  statistically  independent; 

then  it  follows  that  the  interfailure  times  have  a  state-dependent  exponential 
distribution  of  the  form: 

Pr^  >  t}  ■  exp  [-  0(N-i+l)x]  .  (2.1) 

The  resulting  point  process,  a  typical  realization  of  which  is  shown 
in  Figure  1,  is  a  pure-death  continuous-time  Markov  process.  (The  total 
failure  rate  is  a  random  variable  not  only  because  of  the  random  locations 
of  the  failure  epochs  but  also  because  0  and  N  can  be  considered  as 
random  quantities.) 

This  process  can  also  be  viewed  from  a  compe ting-risks  point  of  view, 
in  which  each  error  has  an  independent  probability,  0  dt  ,  of  causing  a 
failure  in  (t,  t  +  dt)  .  The  individual  " lifetimes "  until  an  error  is  un¬ 
covered,  (t^  ;  (i  *  1,  2,  . . .  N) }  are  then  mutually  independent  and  exponen¬ 
tially  distributed  with  parameter  0  .  Our  assumptions  above  then  mean  that 
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all  N  defects  are  considered  to  be  "on  test"  simultaneously  until  failure. 
After  an  interval  (0,t]  ,  we  will  observe  that  a  certain  random  number, 
n(t)  =  n(t)  ,  of  them  will  have  failed  with  observed  lifetimes  {t^,  t2> 

,  and  N  -  n(t)  of  them  will  still  be  "a  failure  waiting  to  occur." 
Since  the  probability  is  Q  =  exp(-0t)  that  a  bug  will  not  surface  in  (0,t)  , 
and  P  =  1-Q  that  it  will,  it  follows  that  the  distribution  of  n(t)  ,  given 
N  and  0  ,  is  Binomial  (N,P): 

Pr{n(t)  =  =  (n  )  (l“e  e  »  (n  **  0,  1»  »••»**)  (2.2) 

with  moments: 

E{n(t)}  =  N(l-e~0t)  ;  V{n(t)}  =  Ne“0t(l-e“0t)  .  (2.3) 

(These  results  also  follow  from  the  solution  of  the  birth-and-death  equations 
of  the  Markov  process  description.) 


P(t) 


h  -H-  t2  t 3  -*■! - 


Figure  1.  Typical  realization  of  overall  failure  rate  over  time 


From  (2.2)  we  can  in  principle  find  the  distribution  of  p(t)  at  any 
t  ;  in  particular,  the  mean  failure  rate  is  the  function 

E{p  (t)  }  =  (N-n(t) )0  =  N0  e"0t  (2.4) 

if  both  0  and  N  are  known ,  and  it  is  this  form  which  suggests  defect 
identification  and  removal  are  very  much  like  reliability  growth  models , 
a  point  already  made  by  early  authors. 

However,  the  software  reliability  problem  is  different  in  the  following 
senses: 

(a)  0  ,  which  reflects  both  the  rate  at  which  certain  portions  of 
the  program  are  exercised  as  well  as  the  probability  that  a 
certain  bug  will  cause  a  failure  to  occur  (plus  possibly  the 
probability  that  an  error  will  be  recognized  by  the  operator), 
is  usually  unknown  a  priori  and  has  to  be  estimated  from  past 
debugging  experience  and  from  the  failure  data  of  this  experiment; 

(b)  N  ,  the  total  number  of  defects  in  the  program,  is  always  unknown 
a  priori; 

(c)  Since  all  failures  will  ultimately  be  found,  the  estimation 
problem  of  primary  interest  after  a  testing  interval  of  t 
hours  is  the  prediction  of  the  distribution  of  unfound  errors, 
nQ(t)  =  N-n(t)  ;  other  possible  measures  of  interest  are  esti¬ 
mators  of  the  local  failure  rate,  p(t)  ,  and  of  the  remaining 
time  until  total  debugging,  t.~^  -  t  . 

By  casting  this  problem  in  a  Bayesian  framework,  similar  to  L-V  and  M-S,  we 
will  see  that  these  distributions  can  be  found  quite  easily.  In  common  with 
M-S,  we  shall  make  the  reasonable  assumption  that  N  is  Poisson  distributed, 
with  parameter  X  hour  *  .  Since  this  only  gives  a  one-parameter  family, 
and  since  it  is  rare  that  X  would  be  known  a  priori,  we  shall  enrich  this 
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assumption  by  additionally  assuming  that  X  is  a  random  quantity  with  a 
Gamma  prior;  this  is  equivalent  to  assuming  that  N  is  a  priori  from  a 
Pascal  (Negative  Binomial)  density,  and  permits  incorporating  large  variances 
in  our  prior  estimates.  In  common  with  L-V  and  M-S,  we  also  assume  that  0 
is  an  independent  random  quantity,  also  Gamma  distributed. 
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3.  FULL-DATA  PRIORS  AND  POSTERIORS 

In  the  usual  case,  the  full  history  of  all  failure  times  will  be 
recorded,  so  that  the  data  from  an  experiment  run  for  t  hours  will  be 
=  t^,  ....  t  ;  n}  where  by  n  we  mean  n(t)  .  For  convenience, 

we  also  switch  to  the  unordered  individual  lifetimes  of  the  n  out  of  N 
failures  that  appear  in  (0,t]  .  It  then  follows  from  first  principles 
that,  given  N  and  0  ,  the  data  likelihood  is: 


P(l>  N,0)  = 


N! 


n 


(N-n) ! 


ir  0e"0ti 


ii-1 


where 


N! 

(N-n) ! 


T  “  T(V  ,N)  =  l  t.  +  (N-n) t 
i=l 


(3.1) 


(3.2) 


is  the  familiar  total-time-on-test  statistic  for  N  items  tested  in  parallel. 
Note  that  (3.1)  remains  the  likelihood  for  any  non-inf ormative  stopping  rule; 
for  example,  instead  of  stopping  after  t  hours,  we  could  also  stop  after 
the  nth  failure,  giving  T=t^  (see  M-S).  Thus,  (n,T)  are  the  usual 
sufficient  statistics  for  0  when  N  is  known. 

However,  in  our  case,  N  is  unknown  and  is  Poisson  (X)  a  priori,  so 

that : 


.N  -X 

p(P.,N  |  X ,0)  =  p(Pt  |  N,Q)  ~f-  • 

Marginalizing  out  the  values  of  N  n  ,  we  get  the  final  data  likelihood, 
given  the  parameters  X  and  0  : 
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p(Pt  i  A,0)  = 


(A0  j  e 


n  -0S  -A  1-e 

P  P 


(3.3) 


where 


S  =  S(P  )  =  I  t  =  l  t  =  l  (n-i+l)T 

i=l  i=i  i=i 


(3.4) 


is  the  total-time-on-test  for  the  discovered  errors.  Thus,  (n,S)  are 
sufficient  for  (A,  0)  . 

It  remains  to  choose  appropriate  priors  for  A  and  0  .  Suppose 
0=0  was,  in  fact,  fixed  and  known  and  0t  was  large ;  then  it  can  be  seen 
from  (3.3)  that  the  Gamma  density: 


p(A)  =  Ga(A  j  a,b)  = 


,a  ,a-l  -bA 
b  A  e 


(3.5) 


would  be  a  convenient  choice  for  the  prior  on  A  ,  since  this  combination 
would  be  closed  under  sampling ,  that  is,  the  posterior-to-data  density, 
P(A  |  V^)  ,  would  also  be  Gamma,  with  revised  parameters: 


a'  =  a+n  ;  b'  =  b+1 


(3.6) 


As  mentioned  previously,  this  is  tantamount  to  assuming  that  N  is  Pascal 
(a,  (b+1)  *)  distributed  a  priori;  the  hyperparameters  a  and  b  can  be 

estimated  from  E{ N }  =  (a/b)  and  ^{N}  =  (a(b+l)/b^)  .  (See  (4.2)  below.) 

-0t 

Conversely,  if  A  =  A  were  fixed  and  all  values  of  Q  *  e  were 
assumed  negligible  compared  to  unity,  then  we  see  that  the  convenient  prior 
for  0  would  also  be  a  Gamma  density,  say  Ga(<t>  |  c,d)  ,  so  that,  posterior- 
to-the-data,  p(0  |  V)  would  still  be  Gamma,  but  with  updated  hyperparameters 


:'  =  c  +  n  ;  d'  =  d  +  S  . 


(3.7) 


The  original  values  of  these  hyperparameters  could  be  estimated  from  the 
relationships  E{0)  =  (c/d)  and  1/ { 0 }  *  (c/d  )  ,  or  from  forming  some  opinion 
about  the  time  for  a  single  error  to  be  found,  which  has  moments  E{0  = 

(d/(c-l))  and  1/ { 0~ 1 }  *  (d2/(c-l)2(c-2))  . 

Now,  unfortunately,  our  likelihood  (3.3)  is  more  complicated,  and  the 

— 0t 

coupling  term,  exp{Xe  }  =  exp{XQ}  ,  means  that,  posterior-to-the-data, 
p(X,  0  |  Pt)  will  have  the  two  parameters  dependent,  even  though  they  were 
a  priori  independent.  Nevertheless,  the  assumption  of  independent  Gamma  priors 
still  turns  out  to  be  the  most  effective  one  for  analytic  simplicity,  and  so, 
using  Bayes'  law,  we  write  the  joint  posterior  on  the  parameters  as: 

p(X,  0  |  Vt)  *  e  XQ  Ga(X  |  a',  b')  Ga(0  |  c\  d')  (3.8) 

where  the  interaction  is  clearly  seen.  Of  course,  as  t  -+•  “  ,  our  posterior 
opinions  about  X  and  0  will  become  independent  again,  since  Q  ■>  0  . 

By  expanding  the  coupling  term  in  powers  of  (XQ)  ,  one  can  write  (3.8) 
in  closed  form  as  the  ratio  of  two  infinite  series,  as  in  Jewell  [4]. 

We  shall  not  pursue  this  here,  but  pass  to  the  more  interesting  problem  of 
predicting  the  unfound  errors. 


PRIOR  AND  POSTERIOR  PREDICTION  OF  UNFOUND  ERRORS 


It  is  well  known  that  the  partition  of  a  Poisson-distributed  random 
integer  N  by  means  of  a  Binomial  process  leads  to  two  independent  Poisson 
processes.  It  follows  that  the  number  of  undetected  errors  at  time  t  , 
nQ(t)  =  N  -  n(t)  is  Poisson  '\Q)  ,  that  is: 


p(nQ  |  A,  0,  t)  = 


fAQ)  °  e-Xt> 


e-0t) 


(4.1) 


If  0  were  fixed,  we  could  use  (3.5)  to  show  that  our  prior  prediction  of 
nQ(t)  would  be  the  marginal  density: 


p(no  I  **‘>  ■  (wq)  {sk)°  • 


(4.2) 


which  is  Negative  Binomial  or  Pascal  (a,Q/(b+Q))  .  Thus,  before  debugging, 
our  prior  opinion  about  the  errors  that  will  remain  after  time  t  would  be 
that  f(nD  I  0»t)  =  (aQ/b)  and  U{nQ  |  0,t)  *  (aQ/b) (l+(0/b))  ,  assuming 
always  that  0  is  known.  In  particular,  as  t  -*•  0  ,  Q  -*•  1  ,  and  we  obtain 
the  Pascal  prior  for  N  previously  mentioned. 

If  we  attempt  to  use  the  Gamma  (c,d)  prior  on  0  with  (4.2)  to  get 
p(nQ  j  t)  ,  the  unconditional  prior  prediction  of  remaining  errors  at  t  , 
we  again  run  into  analytic  difficulties  and  must  settle  for  an  infinite  series 
However,  the  moments  are  no  difficulty,  whence  we  find,  for  example,  the  prior 
expected  number  of  errors  remaining  at  time  t  : 


FUo  I  t)  -  E{»  Q>  -  (f)  (^)C  , 


(4.3) 


and  the  prior  expected  failure  rate  at  time  t  : 


E(p(t)}  -  E{n  0} 

o 


-(f) 


(4.4) 


(d+t) 


M 


.\  •'  ,*  >  '  »  •  1  ■*  *>  *  *  V  W  ^  ^  4~V  ,\  \  \  \  \  V A"  V  V  >, 


Since  p(A,  0  |  P£)  and  P(nc  |  t)  can  only  be  found  in  terms  of 
infinite  series,  it  might  seem  hopeless  to  attempt  to  compute  the  predictive 
distribution  of  nQ  ,  given  the  data.  However,  comparison  of  (3.3)  and  (4.1) 
reveals  that  there  is  a  fortuitous  cancellation  of  the  coupling  term  when 
using 

p(no  I  =//p(n0  I  t)  p(A,  0  |  t>t )  dA  d0  , 

and  we  find  easily  the  predictive  density: 


P(no 


V 


r(a'  +  n  ) 

K  r(a')n  ! 
o 


(4.5) 


where  only  the  normalizing  constant,  K  =  p(0  j  V )  ,  requires  numerical 
computation. 

In  fact,  (4.5)  can  be  simply  computed  by  setting  p(0  |  V ^)  *  1  ,  using 
the  recursion: 


p(nQ  +  1  |  Pt)  = 

P(n0  1  V 


a'  +  nQ\  /  d'  +  tnQ 
b ’  /  \d’  +  t  +  tnQ 


(4.6) 


and  then  renormalizing.  Although  (4.5)  is  not  a  standard  density,  we  see 
from  (4.6)  that,  as  nQ  gets  large,  the  term  involving  c'  and  d'  can  be 
ignored,  so  the  tail  of  p(nQ  |  V  )  is  approximately  Pascal  (a+1,  (b+1)  ^) 
Unfortunately,  in  contrast  to  (4.3),  the  moments  of  (4.5)  must  now  be 


found  numerically. 


In  attempting  to  get  a  point  estimator  from  the  data,  the  usual  idea 

^  «  A 

is  to  find  the  MLEs  (N,0)  from  (3.1).  M-S  have  reported  that  N  can  be 
unstable  for  small  amounts  of  data,  and  might  not  even  exist.  In  our  extended 
model  (3.3)  one  can  find  the  MLEs  (A,0)  from: 


n 

A 


n  .  5  ,  :  .  -0t 
—  =  S  +  A  te 


which  gives  for  0  alone: 


(5.1) 


and  then  A  by  substitution,  and  finally  nQ  *  A  -  n  .  These  estimators 
should  exist  for  all  values  of  t  >  S((?t)/n(t)  ,  but  one  can,  in  fact,  show 

A 

that  (5.1)  gives  positive  0  only  for  t  larger  than  twice  this  critical 
ratio!  In  fact,  solutions  for  the  smallest  permitted  values  of  t  are 
very  unstable  numerically  and  cannot  easily  be  found.  About  the  only  simple 

A 

thing  that  can  be  said  is  that,  for  t  large,  0  approaches  n(t)/S(l?t)  , 
which  is  obvious.  So  much  for  MLEs. 

However,  modes  also  give  useful  point  estimators,  if  one  is  needed. 

We  find  from  the  Pascal  marginal  density  that  the  prior  mode  nQ  =  N  is 
the  smallest  integer  not  less  than  the  solution  to: 


N 


* 


+  1 


a  +  N 
b  +  1 


(5.2) 


From  (4.6),  it  follows  that  the  posterior-to-data  predictive  mode  is  the 
integer  not  less  than  the  solution  to: 
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a'  +  n* 

o 

d'  +  tn* 
o 

b  +  1 

d’  +  t  +  tn* 
OJ 

(5.3) 


which  can  be  found  easily  through  iteration.  Through  arguments  similar  to 

those  in  Jewell  [4],  one  can  show  that  the  first  term  in  (5.3)  is  an  updated 

*  ~  —0  £ 

estimate  N’  of  N  ,  and  the  second  term  is  an  estimate  Q'  of  Q  =  e  , 

so  that  (5.3)  reads  simply  n*  as  N'Q'  .  Limiting  cases  of  the  hyperparameters 

when  the  priors  are  diffuse  or  degenerate  give  various  intuitive  versions  of 

(5.3).  (See  [4].) 

Note  that  (5.3)  exists  for  all  t  >_  0  and  all  data.  However,  for  small 
t  ,  n*  may  be  substantially  less  than  E{nQ  |  1?)  because  of  the  long  tails 
of  the  distribution. 

6.  COMPARISON  WITH  MEINHOLD  AND  SINGPURWALLA'S  MODEL 

Apart  from  the  generalized  stopping  rule,  our  main  difference  with  M-S 

is  that  they  assume  X  is  known,  rather  than  being  a  random  quantity  that 

is  Gamma  (a,b),  a  priori.  The  effect  of  this  can  be  most  easily  seen  in  the 

term  (a'  +  n  )/b'  in  (4.6)  and  the  similar  term  in  (5.3).  The  M-S  model 
o 

can  be  obtained  by  setting  a  *  bX  ,  and  letting  b  -*•  <*>  ;  then  these  terms 

are  replaced  simply  by  X  ,  and  no  longer  depend  upon  the  observations  or 

the  values  of  n  .  At  the  other  extreme,  our  model  can  represent  diffuse 
o 

price  knowledge  about  X  by  letting  b  ■+■  0  (V{X)  -*■  ®)  ;  then  these  terms 

become  simply  nQ  +  n  (or  n* +  n)  ,  and  are  essentially  independent  of  the 

prior.  In  the  general  use  of  the  Gamma  prior,  we  obtain  a  "credibility" 

mixture  between  E{ X }  and  n  +  n  for  this  term. 

o 


’  •  •  .  « •  *  *•..>*>  ,  •*  >  %  VN 

.*  *  V  '  w  •  »  •  *«**v**L*  .'*■  ,  »  . 
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7. 


EXAMPLE 


Suppose  we  set  a  ■  1  and  b  *  0.1  ,  which  corresponds  to: 

E{ A >  =  E{N)  =  10  ;  U{X}  *  100  ;  l/{N)  -  110  ; 
and  then  set  c  =  5  and  d  =  80  ,  that  is: 

E{0}  =  0.0625  V{0}  =  0.000781  E{0_1}  =  20  ;  (/{0-1}  =  133.33  . 

Picking  N  =  20  and  0  =  0.025  ,  one  realization  of  ordered  failure  times 


(k  *  1,  2,  . . 

•,  20)] 

gave : 

3.282 

3.889 

5.900 

5.903 

5.996 

6.055 

9.829 

10.40 

10.95 

11.90 

11.96 

14.40 

19.98 

20.23 

20.38 

22.91 

30.05 

35.35 

35.69 

59.44 

*  *  * 

Figure  2  shows  the  integral  solutions  n  and  N  *  n(t)  +  n  to  the 

o  o 

*  * 

predictive  mode  equation  (5.3),  versus  t  .  For  t  <  t^  ,  both  nQ  and  N 
are  zero  because  no  failures  have  occurred,  even  though  our  prior  expectation 
is  F{N}  =  10  .  Thereafter,  the  curve  jumps  upward  at  each  failure  epoch, 
followed  by  integer  decreases  as  t  increases  without  further  events  (Figure  2 
shows  only  a  straight-line  approximation  to  the  true  behavior) . 

Figure  3  shows  the  exact  predictive  means  E{nQ  I  ^t)  and  F{N  |  V 1  } 
for  t  *  0(5)80;  these  were  computed  numerically  from  the  full  distributions. 
In  general,  these  behave  more  smoothly  than  the  posterior  modes  and  seem  to 
converge  towards  the  true  N  for  lower  values  of  t  .  Further,  for  small  t 
they  are  influenced  by  the  prior,  rather  than  starting  at  zero.  (Figure  3 
shows  only  a  straight-line  approximation  to  the  true  behavior.) 

However,  the  main  result  is  not  these  point  estimators,  but  the  fact  that 
the  full  predictive  density  can  be  obtained  from  (4.6).  Figure  2  shows  the 
(Pascal)  density  obtained  at  the  start  of  the  experiment  (t  ■  0);  for  our 


.V 


parameters  the  initial  mode  is  zero.  After  a  couple  of  errors  are  found,  the 
mode  shifts  to  a  positive  number  and  the  variance  begins  to  reduce,  giving 
curves  similar  to  Figure  5,  where  t  =  10  and  n(10)  =  7  .  Thereafter,  the 
curve  essentially  remains  unimodal,  with  nodes  and  means  given  by  Figures  2 
and  3,  with  the  variance  decreasing  almost  continuously.  After  about  t  =  46 
n(46)  =  19  ,  the  predictive  density  again  has  its  mode  at  the  origin,  with 
mass  on  only  a  few  values  of  nQ  ;  Figure  6  shows  a  typical  density,  at  t  =  7 
10  time  units  after  the  last  error  has  actually  been  found. 

In  this  way,  we  see  how  easily  one  obtains  full  predictive  information 
about  undetected  errors  for  use  in  decision  problems  or  in  making  performance 
guarantees . 


Figure  2.  Approximate  behavior  of  predictive  inodes,  n  and  N  ,  versus  testing  interval 
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